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A COMBINED NEWTON- RAPHSON AND GRADIENT 


PARAMETER CORRECTION TECHNIQUE FOR SOLUTION 
OF OPTIMAL- CONTROL PROBLEMS* 

By Ernest S. Armstrong 
Langley Research Center 

SUMMARY 

A parameter correction technique is developed to solve a boundary-value problem 
which frequently occurs in optimal- control theory. It is assumed that an indirect 
optimal- control method has been applied to a controllable dynamic system with a two- 
point boundary-value problem resulting such that the boundary conditions take the form 
of a set of unknown parameters to be determined to meet an equal number of terminal 
conditions. The optimal- control law is a piecewise continuous function with discontinu- 
ities occurring only at the zeros of certain continuous functions. A procedure is devel- 
oped to improve upon an assumed set of parameters so that, by repetitive use of a cor- 
rection formula, a monotonic decreasing sequence of values of a positive definite function 
that measures the terminal errors is produced. The direction of the correction vector 
is found to lie between the directions given by the gradient and the Newton- Raphson 
procedures. 

Integral equations are derived for influence matrices that describe the effect of a 
change in the parameters on the terminal conditions. 

The procedure is successfully applied to the determination of both planar and non- 
planar fuel-optimal trajectories for a space vehicle which is launched from the surface 
of the moon and required to rendezvous with a space vehicle in a circular orbit. 

INTRODUCTION 

In recent years, control theory has been expanded to include the area of system 
optimization. This expansion has brought about a new design philosophy. Control 

*This report is based in part upon a thesis entitled "An Algorithm for the Iterative 
Solution of a Class of Two- Point Boundary- Value Problems Occurring in Optimal- Control 
Theory" offered in partial fulfillment of the requirements for the degree of Doctor of 
Philosophy in Mathematics, North Carolina State University of Raleigh, Raleigh, North 
Carolina, June 1967. 



functions may now be chosen to optimize in some given sense the system response to the 
control action; for example, a control law may be found to force a system to a given state 
while some functional of the system variables is minimized. This new area of research 
is termed optimal- control theory. 

Present-day methods of calculating optimal- control solutions can be grouped into 
two classes: direct and indirect. Both methods are designed to minimize the value of 
some functional. A direct method depends upon a comparison of the values of the func- 
tional at two or more points. An indirect method is used to find a solution by means of 
necessary (and sometimes sufficient) conditions for a minimum. Typical direct methods 
are contained in references 1 to 3. Necessary conditions to be used in an indirect 
approach are found by applying the Pontryagin maximum principle (ref. 4), or the calculus 
of variations (ref. 5), or dynamic programing (ref. 6). In general, the necessary condi- 
tions take the form of a set of nonlinear differential equations with both initial and final 
boundary conditions; that is, in order to obtain explicit solutions, a nonlinear two-point 
boundary-value problem must be solved. 

The advent of high-speed computers has made feasible the solution of optimization 
problems by the method of successive approximations. This procedure is markedly 
illustrated by the success of the aforementioned direct methods. In these methods, a 
control history is first assumed and then successively improved upon by the computation 
of time -dependent corrections arrived at through the use of gradient (refs. 1 and 2) or 
conjugate-gradient (ref. 3) theory in function space. Although many useful results have 
been obtained in this manner, direct methods, in general, do not guarantee that the solu- 
tions obtained satisfy the necessary conditions of the indirect theory. 

A more rigorous, but computationally more difficult, approach is the use of neces- 
sary conditions of the indirect theory for the generation of optimal results. In this way, 
one of the theories of references 4, 5, or 6 is applied, and then an attempt is made to solve 
whatever boundary-value problem might ensue. This approach is adopted herein. 

The purpose of this report is to present a successive approximation procedure for 
attacking a class of two-point boundary-value problems which frequently occurs in the 
application of indirect optimization theory. Basically, the boundary-value problem is 
one in which the optimal- control law is piecewise continuous and in which there are a 
number of system parameters to be determined to meet an equal number of terminal con- 
ditions. A parameter correction procedure is developed in which an assumed set of 
parameters can be improved upon so that, by repetitive use of a correction formula, a 
monotonic decreasing sequence of values of a positive definite function that measures the 
terminal errors is produced. The direction of the parameter correction vector lies 
between the direction given by the gradient and the Newton- Raphson procedures (ref. 7). 
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Integral equations are derived, the solutions of which yield influence matrices that 
describe the effect of a change in the parameters on the terminal conditions. 

In order to exemplify the usefulness of the procedure, the Pontryagin maximum 
principle is applied to determine planar and nonplanar fuel- optimal trajectories for a 
space vehicle which is launched from the surface of the moon and required to rendezvous 
with a space station in a circular orbit. The technique is then successfully applied to 
solve the resulting two-point boundary -value problem. 

SYMBOLS 

A,D,M,N,K,L,S,G constant matrices 

bi positive weighing elements (i = 1, 2, . . . m) 

B m-dimensional diagonal matrix with elements bj 

\Ib m-dimensional diagonal matrix with elements \jb[ 

c effective exhaust velocity 

Ci elements of c (i = 1, 2, . . . m) 

C used as a(t)|£ to designate continuous part of function a(t) 

d = xp 2 (x 1 + R Sx ) + i// 4 (x 3 + R Sy j + «//g (x 5 + R Sz ) 
e m-dimensional vector with elements ei 


e i 

*(*'*t) 

f 

f i 

f Q (x,u) 


terminal errors (i = 1, 2, . . . m) 

m , o 
e • Be _ V ^i e i 

2 L 2 

i=l 

n-dimensional column vector with elements fj 
function introduced in equation (la) (i = 1, 2, . . 
integrand of F 


n) 
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function to be minimized, f Q (x,u)dt 

*•0 

n-dimensional column vector with elements defined in equation (2) 


g^ elements of vector g (i = 1, 2, . 

h = co\^P v \p 2 , . . . 
n 

S = X ^k*k 


• n) 


k=0 


H(a) = |(1 + sgn a) 
i,j,k unit vectors 

I identity matrix 

J(t) set difference, |^ 0 ,tj “ S(t*) 



l 

m 

m(t) 

m 0 

max 

M 

n 

P 


total number of switching functions 

number of unknown parameters and terminal conditions 

total vehicle mass 

initial mass of launch vehicle 

maximum 

maximum value of H with respect to choice of u 
dimension of x and \p in equation (2) 

smallest positive integer where p^(a, t) i= 0 
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r 


dimension of u 


I* — Ry — Rg 

r x ,r y ,r z elements of ? 

r x ,r y ,r z elements of r 

_ \ 1/2 

R s satellite orbital radius about moon, ( R s • R S J 

Rs = i R s x + l^Rsy + kRs z 
Rs x » R s y > R s z elements of R s 

Rg x ,R Sy ,R Sz elements of R s 

— o—o — — 

R s ,R S initial values of R s and R s , respectively 

Ry magnitude of position vector of interceptor vehicle, (Ry • Ry) ^ 

R v = i Ry x + f Ry y + kRy z 
R Vx ,R v y ’ Rv z elements of Ry 
Ry x ,Ry y ,Ry z elements of Ry 

— O -L. O — 

Ry ,R V initial values of R v and R v , respectively 

s dummy integration variable 

r 1 (a > 0) 

sgn a general signum function defined by <-l (a < 01 

^Unspecified (a = 0) 

sgn Pj a particular signum function defined in equation (4) 

S(t*) set of switching points of all switching functions p (q = 1, 2, . . . l) 

4 

S n (t*) set of switching points of switching function p 

t element of ("t 0 ,t^| 
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tf final time 

t Q initial time 

t + ,t" t approached through values larger than t and smaller than t, respectively 

t* arbitrary switching point 

t] 1 ' ith switching point in S(t*) 

T magnitude of thrust- control vector 

T thrust- control vector 

Tl constant matrix defined in equation (B4) 

T 2 (t) matrix defined in equation (B5) 

u r- dimensional vector with elements uj 

u unit vector using some elements of u 

Uj control elements (i = 1, 2, . . . r) 

U r-dimensional Euclidean space containing u 

v total number of switching points in S(t*) 

v = COl^Xj, x 2 , . . . x 6 J 
v 0 initial value of v 

V n- dimensional Euclidean space containing x 

x,y,z coordinates of axis system in figure B-l 

^ + R s x ) 2 + (x 3 + Rs y ) 2 + (*5 + R s z ) 2 

x vector with elements xj 
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state variables given by equation (la) (i = 1, 2, 


c o = f oM) ds 


x r ,y',z' coordinates of axis system in figure B-3 


x° = x(t Q ) 


x-*- = x(tf) 

X,Y,Z coordinates of axis system in figure B-3 
x(a,tf^ equation (15) evaluated at (a,tf) 


vector defined in equation (Bll) 


m-dimensional vector with elements 


unknown parameters (i = 1, 2, . . . m) 


vector defined in theorem 1 (i = 0, 1, . . .) 


variable converting 6a • 5a ^ iP into 5a • 5a - v z + fi* = 0 


ts\i ,.2 , q2 


largest allowable value of thrust magnitude 


6a, Aa increment in a 


increment in a in gradient direction of -E^a,tf) 
5 Nr increment in a in Newton-Raphson direction 


6v = A 6 a 


AE(a,tf) = E^a + 6a, tf) - E^a,t^ 
AE^a,tj^ defined by equation (7) 


arbitrary m-dimensional vectors 


t 



angles defined in figure B-2 


8 C ,<P C 

Ov°,(Py° angles defined in figure B-4 
i 0 ,6 0 ,cp 0 angles defined in figure B- 3 


\(u) 


Lagrange multiplier 


value of X associated with v 




0 p Oa 

ith eigenvalue of matrix ~^=r B — (i = 1, 2, . . . m) 
gravitational parameter of moon 


bound on magnitude of 6a 


bound on magnitude of 6a 1 


P 

Pi 


(p 




l -dimensional column vector of elements 
switching function (i = 1, 2, . . . l) 
arbitrary value of te [to’tf] 

scalar function used as stopping condition 


\|^= ]fp 2 2 + ^ 4 2 + 2 


^i 

^(a,tf) 


u> 


a> = ka> 




vector with elements i/a (i = 1, 2, . . . n) 

variables introduced by Pontryagin maximum principle (i = 0, 1, 2, 
equation (17) evaluated at 

angular velocity of moon about axis of rotation 
angular velocity of target in orbital plane 


• n) 


absolute value of a 
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I 


I I 


i 


a designates a is vector 

|| a|| = (a • a) 1 / 2 


[a,b] 

(a,b) 

a(t) 

a(t) 


closed interval 
open interval 

first derivative of a(t) with respect to t 
second derivative of a(t) with respect to t 


a^(t), ^ a ft) nth derivative of a(t) with respect to t 


dt 1 


n 

b = aibi if a = col^a^, . . . a n ^ and b = col^bj, 


i=l 


l-l 


P q(t*) 

9x 

8 y 


9 2 a 
9y 9x 


sum over all switching functions p (q = 1, 2, 

, 4 

switching points at t* 


9X: 


bn) 

. 1) which have 


M x N Jacobian matrix with elements -5— i where i = 1, 2, . . . M and 


9 Yi 


j = 1, 2, . . . N if x = col^Xj, 


x m) and y = coi(y r • • • y N ) 


M X N matrix with elements 


9 2 a 


9xi 9yj 


7 where i = 1, 2, . . . M and 


j = 1, 2, . . . N if x = col^Xj, • • • X M ) and y = col(y x , . . . y N ) 
null vector 


e belongs to a set 

= defined by 

~ over a variable indicates a function of a and t obtained by substitution for 

x(a,t) and i//(a,t) into a function of x(a,t), if/(a,t), a, and t 

~ over a variable indicates a function of t obtained by substitution for x(a,t) 

4/(a,t), and a in a function of x(5,t), a, and t 

' over a matrix denotes matrix transpose 
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Subscripts: 


j,k,m,n,q,r,v,M,N integers 

Superscripts: 

i,p integers 


PROBLEM STATEMENT 
Indirect Optimal- Control Theory 

Consider the behavior of a dynamical system the state of which at any instant of 
time is characterized by n variables Xj, Xg, . . . x n . For example, these variables 
might represent the position and velocity coordinates of a space vehicle. The behavior of 
the system is simply the time variation of the vector variable x = col^Xp Xg, . . . x n ^, 
commonly referred to as the state vector. The vector space V of the vector variable x 
is termed the state space of the system. 

Assume that the state of the system can be controlled; that is, a set of system inputs, 
the manipulation of which governs the state, are available. Assume that there are r 
such controls and that they are characterized by a point u in a region U of 
r-dimensional Euclidean space. For the purposes of this report, 

u(t) = colju^t), . . . u r (t)] is said to be admissible if each component is a piecewise con- 
tinuous function such that u(t) belongs to U for each t (t 0 2 t ^ t^. The initial time 
t 0 is assumed fixed, but the final time tf may be either free or fixed. Let the behavior 
of the dynamical system be characterized by the autonomous differential equations 

dxj . , 

= fi( x i> x 2 > ' * * x n> u 2 > • • • u r) (i = 1, 2, . , . n) 

In vector form 

= f(x,u) 

where f(x,u) is a vector function of x and u. The functions fp for every xeV and 
ueU are assumed to be continuous with respect to all variables x^, . . . x n and 
Up . . . u r . Also the functions fj are continuously differentiable with respect to 
Xp . . . x n ; that is, 

8f i 

f i (x, u) and (i = 1, 2, . . . n; j = 1, 2, . . . n) 

are defined and continuous for all xeV and ueU. 
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The admissible control u(t) is said to transfer the point x° to the point x* if 
the solution x(t) of ^ = f (x,u) (x(t Q ) = x°) passes through the point x* at a certain 
instant of time tf ; that is, x(tf) = x 1 * . 

The quality of system performance is now assumed to be measured in terms of 
some functional 

F = J ^ f 0 [x(t),u(t)]dt 
*o 

For example, the smaller the value of F, the better the system behaves. The scalar 
function f 0 [lc(t),u(t)] is defined and differentiable with respect to Xf (i = 1, 2, . . . n). 

The optimization problem to be considered is the following: Given the dynamical 

system ^=f(x,u) with x(t Q ) = x° and a point x*eV, find an admissible control ueU 

(if any exist) which transfers x° to x* such that F takes on the least value. In gen- 
eral, necessary conditions for the solution of this problem are given by the Pontryagin 
maximum principle (ref. 4). 

In order to formulate the maximum principle, define x Q (t) such that x Q (t) = f Q (x,u) 
(x 0 (t 0 ) = o). (Note that x Q (tf) = F.) In addition to the system 

dxj , „ 

— = fj(x,u) (i = 0, 1, 2, . . . n) (la) 

consider another system of equations 

dj£/, y' 9 fi 

"dT = " Z (i = 0, 1, 2, . . . n) (lb) 

3=0 1 

in the auxiliary variables \p Q , i//j, . . . xf/ n and define 

n 

S = I 

i=o 

For fixed values of xj and (i = 0, 1, . . . n), H becomes a function of ueU. 

The Pont ryagin maximum principle .- Let u(t) (t 0 5 1 1 tf) be an admissible con- 
trol which transfers x° to x* by equation (la) such that x Q (tf) is minimized. For 
this case, it is necessary that there exist a nonzero continuous vector 
[*o(t), ^f(t), . . . i^ n (t)J satisfying equation (lb) such that: 

(1) The control u(t)eU maximizes H for fixed x A and (i = 0, 1, . . . n) 

at the point u(t) ; that is 
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~ n 

[fx i (t),i// i (t),u ; j(t)| = max H = max 
J ueU ueU 

I V# ”) 

_k=0 




(i = 0, 1, . . . n; j = 1, 2, . , . r) 

Then, M^f(t),Xj(t)J represents the maximum value of H attained by substituting 
u = u(t) into H. 

(2) For all te|t 0 ,tjj, ^ 0 (t) = Constant = 0 and M|xf(t),i//^(t)| = Constant = 0. 

The statement of the Pontryagin maximum principle is for an autonomous dynamic 
system with x° and x* given and tf undetermined. Extensions of this theorem to 
autonomous systems with tf fixed and to nonautonomous systems with tf free and 
fixed are given in reference 4. Cases in which some x° and x* are unknown involve 
another facet of the theory known as the transversality condition (ref. 4). For illustra- 
tive purposes, consider the following example: 

EXAMPLE: 

Let the dynamical system be characterized by 

AS + DS 
dt 

where x is an n-dimensional column vector, u is an r-dimensional column vector, 
and A and D are (n X n) and (n x r) matrices, respectively. Constrain the controls 
so that 

|uf| il (i = 1, 2, . . . r) 

Find Uf(t) such that the dynamical system is transferred from x° to x 1 in minimal 
time; that is, x Q (tf) = tf or f 0 (x,u) = 1. 

From the Pontryagin maximum principle, the optimal control u(t) maximizes 

H=^ 0 + ^'5f=^o + *** + 

where if/ Q = 0, \p is the n-dimensional column vector i// = col^j, . . . tyn), and the 
prime denotes transpose. 

Obviously, H is maximized for u(t) = sgn^xf/ d] with t 0 £ t S tf where 


The system 


r 1 

sgn a = / -1 

^Undefined 

~(t) = Ax(t) + D sgn|j^'D) 


(a >0) 
(a < 0) 
(a = 0) 
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^(t) = -A>(t) 

with the boundary conditions x(t 0 ) = x° and x(tf) = x* now results. From condition (2) 
of the Pontryagin maximum principle 

i£'(t)jAx(t) + D sgn[^'(t)D]U + = 0 

Hence, there exist (2n + 1) conditions for determining the variables x(t), \p(t), and tf. 

Note that the form of the optimal control u(t) follows readily from the maximiza- 
tion of H. This feature is desirable. However, the optimal control is governed by i//(t) 
and the initial conditions ^(t Q ) are not given. Thus, there exists a nonlinear two-point 
boundary -value problem which can be stated in the following form: Determine the (n + 2) 
unknown parameters \p 0 = 0, i^(to), and tf such that at tf, the (n + 1) terminal condi- 
tions x(tf) = x* and Mjx(tf),^(tf)J = 0 are met where x(t) and \p(t) satisfy 

gf(t) = Ax(t) + D sgn(y/Dj (x(t 0 ) = x°) 

ff(t) = -A'#) 

Because must be a constant greater than or equal to zero, the boundary -value prob- 

lem can be separated into two cases. Both cases involve (n + 1) unknown parameters to 
be found so as to meet (n + 1) boundary conditions. In one case, is set equal to zero; 

in the other case, i p Q is chosen as some negative constant. 

Such boundary-value problems typically result from the maximum principle and 
other indirect theories and are characteristic of their main difficulties. Generally, 
because x Q (t) is completely specified by 

x 0 = f 0 (x,u) (x 0 (to) = 0) 

and x(t) and u(t) are determined, x Q can be eliminated from the boundary -value 
problem. The corresponding auxiliary variable i// Q can be eliminated by separating the 
problem into two cases as in the foregoing example. Thus, any two-point boundary -value 
problem originating from the maximum principle can be made to involve only the differ- 
ential equations for Xj and i/q (i = 1, 2, . . . n). 

A Particular Boundary-Value Problem 

The purpose of this report is to present a successive approximation procedure for 
attacking a class of two-point boundary-value problems which frequently occurs in indirect 
optimization theory. The particular class of boundary-value problems to be considered 
and the mathematical assumptions concerning it are now presented. 
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The general result of applying an indirect method such as the Pontryagin maximum 
principle is a set of necessary conditions which can be arranged as 2n differential 
equations with mixed-boundary conditions. With the differential system defined over 
|t 0 ,t^], the boundary-value problem to be considered is one in which t Q is known, the 
optimal control u(t) is piecewise continuous, and there are (m ^ 2n) parameters repre- 
sented by the column vector 3 = (a^, . . . o; m y to be chosen such that m terminal 
conditions are met. The parameters are some or all of the initial values of the differen- 
tial equations and possibly the duration tf - t 0 . By writing x and \p as x(3,t) and 
i//(3,t), respectively, to indicate their dependence on 3, the terminal conditions can be 
represented as 

e i[S(a,tf),i/7(3,tf),3,tfJ = 0 (is 1, 2, ... m) 

The 2n differential equations can be written in the form 
^x(3,t) = f[x(3,t),u,3,t] 

J^(5,t) = -(||) jj(3,t),u,5,t]i£(3,t) = §[x(3,t),^(3,t),u,3, 



Assume that the optimal-control functions take the form 

u = ujic(3,t),i/7(3,t),sgn p(i(3,t),^(3,t),3,tj,3,t| (3) 

where sgn p is an l -dimensional column vector with element sgn p^. The p^(t) 

(i = 1, 2, . . . 1) are continuous functions of t for given x(3,t), i//(3,t), and 3. Then, 
sgn p^(t) is defined as 



The functions u, f , and g are to be continuous in x, t p, u, sgn p, and 3 and 
piecewise continuous in t with points of discontinuity occurring at those values of t 
for which p^(t) = 0 (i = 1, 2, ... l). 

Assume that: 

(1) For given x(t), i^(t), and 3, the zeros of p^t) are finite in number where 
i = 1, 2 , 1. 
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(2) These zeros vary continuously with a. 

(3) If t* is a zero of p^(5, t) forgiven x(a,t) and \p(a,t), then p^(a, t) is 

assumed to be continuously differentiable in t at t* of order equal to the first non- 
vanishing derivative from the left of p^(cf,t) at t*; that is, if p is the smallest posi- 
tive integer such that p^(a,t*-) =* 0, then p^(ci,t*“) = p^(a,t*). Also, in such a case, 
-g^p-Pj(5, t) (i = 1, 2 , 1; j = 1, 2, . . . m) is assumed continuously differentiable in 


c ] 

t at t* of the order p - 1. 

Mi 9gi 

(4) The matrices -5 — and - — 
8oi k 8oi k 


(i = 1, 2, . . . n; k = 1, 2, . . . m) are bounded 


and continuous with finite interior limits as the boundaries are approached on the set 
J(tf) = [t 0 ,tfj - [The set of switching points of pj(t) (j = 1, 2, . . . Z)J. Switching points 
are particular zeros of pj(t) as discussed subsequently. 

An iterative procedure is now developed to improve upon an assumed value of 5 
in such a way that, when equations (2) with u given by equation (3) is satisfied, ej = 0 
(i = 1, . . . m) results. 


THEORETICAL DEVELOPMENT OF CORRECTION TECHNIQUE 


Iterative Logic 

Given a value of a, equation (2) can be integrated and a set of values of 
ej(x(5,tf),i/'(5,tf),3,tf^ (i = 1, 2, . . . m) can be computed. In order to simplify the nota- 
tion, let these values of ej be represented by e^a,tf^. 

Let e(a,tf) be a column vector with elements ej(a,tf) (i = 1, 2, . . . m). Define 
a function E^5,tf) as 

E(a,tf) = - 

where B is an (m x m) positive definite diagonal matrix. The scalar function E(a,tf) 
can then be used as a measure of closeness to a solution because finding an a for which 
E(c5,tf) = 0 is equivalent to finding an a for which ej(3,tf^ = 0 (i = 1, 2, . . . m). 

The elements of B are to be used as weighing coefficients. 

Now assume a value of <3; for example, 5°. Integrate equations (2) and evaluate 
E(a°,tf). If E(a°,tf) = 0 for all practical purposes, then the boundary -value problem 
is solved; if not, let a° be changed by an amount 5a°. This change causes e(a°,tf) 
to become e(a° + 6a°,tf). Assume that 


(“,tf) ‘ Bi(a,tf) 




15 


Therefore 


“=■ °» A (a! °, tf) , l// (a °, tf) , a , tfj^ 


exists for row i = 1, 2, . . . m and column j = 1, 2, . . . m. 

A -/-< 


Ae^a °,tf] = e(a° + 6a 0 , tf) - e(a °,tf) = ^£(a°,tf)6a° + o(6a°) 


( 6 ) 


where 


11 °( 6q!0 ) 


lim il - 
6a °| 1—0 ll 6a! 


7T. O 


= 0 


The change Ae(a°,tf) causes a change AE(a°,tf) in E(a°,tf) as given by 
AE(a°,tf) = E(a° + 6a°,tf) - E(a°,tf) = Be(a°,tf) • Ai(a°,tf) + i Ae(a°,tf) • B Ae(a°,tf) 

Let AE(a°,tf) represent the value of AE(a°,tf) obtained by replacing Ae(a°,tf) by 
the first-order term in 6a ° of equation (6) so that 

AE(a°,tf) = ■^•(a°,t f )Bi(a°,tf) • 6a° + | 6a° • ^-(a°,tf)B ||(a°,tf)6a° (7) 


It will be assumed that ||6a°|J can be made so small that the behavior of AE(a°,tf) 
can be gathered from AE(o°,tf). In particular, ||6a°|| is so small that the algebraic 
sign of AE(a°,tf) is the same as that of AE(a°,tf). Analytically, this smallness on 

||5a°|| can be represented by ||5a°|| = (v>0). 

In appendix A, lemmas 1 and 2 establish that unique solutions 5a°(v) and 
\(v) <0 of the system 

,2 


,«o||. = „2 


6a 


6a° = - 


% W - “J ' 


1-1 

°e(-.o i.a \ T I 9 e fco 




30 ! 


( 8 ) 


exist if v is sufficiently small (lemma 1). The solutions maximize the absolute value 
of AE(a°,tf) subject to the conditions ||5a°|| g and AE(a°,tf) < 0 (lemma 2). 

Also, £E(a°,tf) is negative definite in 6a° for 6a° to satisfy equations (8) (lemma 3). 
Thus, given a 0 and v where AE^a°,tf) approximates AE(a°,tf) and lemma 1 is 
satisfied, the replacement of a 0 by a' = a° + 6a°, where 6a° satisfies equations (8), 
yields E(a',tf) < E^a°,tf). This property follows from lemma 2. Repetitive use of this 
procedure generates a monotonic decreasing sequence of values of E(a x ,tf) with 
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a* + * = a* + 8a* (i = 0, 1, . . .) and 6a* satisfying equations (8). The values of v for 
each i may not be the same. Because E(a*,tf) = 0 (i = -0, 1, . . .), the sequence is 
bounded from below and therefore converges. The point of convergence is where 
SE(a*,tf) vanishes which, by lemma 3, occurs where 5a* = 0. Thus, the following itera- 
tive logic is established. 


Theorem 1 .- Given values of a and v, for example, a u and v Q , respectively, 
for which AE(a°,tf) approximates AE(a°,tf) and lemma 1 holds, then the use of 

(i = 0, 1, . . .) with 6 a*(jq) and |\(iq)| given by the simultaneous 


ai+1 = q4 + 5a;i 
solution of 


| 63 i || 3 = n 


_ „ 2 


( 9 ) 


and 


6a* = 


^( ai ^) B H^^f) + | X (^i)| I J ^(a i ,tf)Be(a i ,tf) (10) 


3e 


1-1 


generates a monotonic decreasing sequence E(a*,tf), if at each point of the sequence, 
values of iq can be found for which XE(a*,tf) approximates AE(a*,tf) and lemma 1 
• holds. The sequence E(a*,tf) converges to a value of E^a,tf) for which 6a of equa- 
tion (10) vanishes. 


The determination of values of iq for theorem 1 can be simplified by manipulating 
|\(iq)| in equation (10) until AE(a*,tf) approximates AE(a*,tf) and then by computing 
iq from equation (9). A better approach is obtained by manipulating the |\(iq)| inequa- 
tion (10) until 

AE(a*,tf) = E(a* + 6a*, tf) - E(a*,tf) < 0 


because theorem 1 may only establish sufficient (and not necessary) conditions that equa- 
tions (9) and (10) generate a monotonic decreasing sequence E(a*,tf). 


Convergence 

From equation (10), a limit point is reached when 

i, tf) B e (a *, t^ = 0 


9a 


The existence of the inverse of -^r(a*,tf) determines whether e(a*,tf) = 0. Obviously, 
the existence of ^!(a*,tf) * is sufficient but not necessary for E(a*,tf) to vanish when 

_ T _ T 

|lr(a*,tf)Be(a*,tf) does. The vector -rp-(«*,tf)Be(a,tf) can be recognized as -^-(a,tf) 


so that the arrival at -^p-(a,tf)Be(a,tf) = 0 means that an extremal of E(a,tf) in the 
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hyperspace a,E(5,tf) has been reached. The possibility of e(a!,tf) #0 at such a point 
always exists unless E(a,tf) is strictly convex in a; that is, the Hessian matrix 

02-pi 

— =-(o? ,tf ) is positive definite. Because 
da 2 ' ' 


E(3 , tf ) . i 

i=l 


then 


m r. 


0(a,t f ) = l + eiKt f )^ a ,t f ) 


8^E — 

In general, it is not analytically evident that ^«{a,tf) need be positive definite. 

da 2 ' 

Also, the computations necessary to evaluate 3^E/9 numerically and to test for posi- 
tive definiteness would be excessive. It is therefore recommended that the procedure be 
used without attempting to examine the definiteness condition. However, note that when 
||e(«, tf )|| is so small that 

^ . v* 3ei'/_ , 9 ei 


o m , . 


i=l 


then 


9 m r -i2 


when 4 is an arbitrary column vector; that is, E(a,tfJ is locally convex when 
||e(a,tf)|| is small. 

A procedure which relies on the numerical inversion of even a theoretically non- 
singular matrix may incur stability problems. The matrix may be numerically near 
singularity in the sense that a substantial loss of significant figures results from the 
inversion process. The algorithm has the advantage that the nonsingularity of 

■^=r(a,tf)B !=(a , tf^ + | x|l can be controlled by manipulation of | x|. Even though 

•jp(<5,tf) is singular, -^r(a,tf)B |^(a,tf) + j Xjl can be made strongly nonsingular by 
increasing J X |; thus numerical stability is provided by the iteration process. 


18 



Relation to Gradient and Newton- Raphson Processes 


An important feature of the algorithm is observed in the limit as I A I either 
increases or decreases. As A 


increases. 


and from equation (10) 




da 


approaches 




I, 


= - \t\ = 6 “c 


(ID 


Equation (11) is simply the well-known method of gradients (ref. 7); that is, 6a is 
chosen to lie in the direction of the negative gradient of E(a,tf) with respect to a at 

(«’*)• K 1 exists, then as |a| decreases 


6a — - 



NR 


which can be recognized as a Newton- Raphson iteration procedure. For arbitrary A, the 
process designated by equation (10) represents a combination of the gradient and Newton- 
Raphson techniques. 


The result given by equation (10) is similar to that obtained in reference 8. Ref- 
erence 8 shows that the angle between the vector 5a of equation (10) and the gradient 
vector of equation (11) is a continuous monotonic decreasing function of | A| such that 
as | A | — °°, the angle approaches zero. Because the gradient vector of equation (11) is 
independent of | A j, it follows that 5a of equation (10) rotates toward the 5a of equa- 
tion (11) as f A| - °°. The direction of the correction vector of equation (10) then lies 
between the directions given by the gradient and Newton- Raphson procedures; that is, 6a 
given by equation (10) is between 6a^ and 6a in the sense that 6a always makes 
a smaller angle with the gradient direction of equation (11) than does 5a 


The Matrix -^=(a,tf) 

9a' 

q" 

The form of the matrix -^(a,tf) depends upon the nature of the final time tf in 
the boundary-value problem. These forms can be divided into the following cases: 

Case 1, where tf is a constant.- Because e|^(a,tf),i//(a,tf),a,tfj may depend 

implicitly, as well as explicitly, upon a through the presence of x(a,tf) and 1//(a,tf), 
then 
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^0 ’ r~ ^ — i 



+ 



^0-j f- V 


(j = 1, 2, . . . m; k = 1, 2, . . . m) (12) 

Case 2, where tf is unspecified.- If the final time tf is unspecified, it can be 
treated as a quantity which can be initially guessed and then corrected to meet the termi- 
nal conditions. 


When the final time is unspecified, an (m - 1) dimensional parameter vector a 
and tf must be determined to meet m terminal conditions. Adding an mth column 

■gjffp, tf) , P («. tf) , 5 , tf) to the (m x m - 1) matrix |I(a,tf) yields an (m x m) matrix 

which, when used in equation (10), generates a correction vector where the first (m - 1) 
elements represent 6a and the mth element represents 6tf. By this method, a correc- 
tion to tf can be computed which is influenced by the other 6a. In the newly formed 
(m x m) matrix, the first (m - 1) columns have elements given by equation (12). The last 
column is given by 


dt f 


86ir "(a,tf), xl/(a, tf), a, tf] = ^j^(a,t),^(o!,t),a!,t| 


t=tf 


+ 



9e 




0 - 1, 2, ... m) (13) 


The last two terms of equation (13) regard the final time as being fixed and take 


c(of,tf^ and i//(<r,tf) as a result of a change in tf at t Q . 
Thus, ^(“.tf) and , tf) satisfy the same type of equations as do and 


into account any changes in 

9 X/-x\ 


da i. 


(a , t() in equation (12). 


Case 3, where tf is implicitly specified.- A situation may occur in which the final 
time tf is free to vary from iteration to iteration but is determined after a choice of a 
is made through 



This case may be treated in one of the following manners once a tf where 
<p|lc(a , tf) , ip(a , tf) , a , tf = 0 has occurred: 

(a) Choose in equation (10) so large that the change in tf from a to 

o' + 6o! is small. Then, use the method in case 1 but iterate at the final time deter- 
mined by (p . 

(b) Adjoin ,tf),«,tfj to e£x(<5,tf), 4/(a, tf), a, t^j and tf to a and 

use the method as in case 2. 

(c) Add to equation (12) the term 



dfCt) = ^ <p(x(a,t), !//(«, t),a,t| 


t=t f 


* 0 


9x- 8;//. 

and take into account (<5,tf) in the computation of and . 


Cases 1 to 3 characterize — (a,tf). In order to complete the representations for 


|=(a,tf), it is necessary to compute |=(<*,tf) and tf) for t f treated as if it 


were fixed (cases 1, 2, 3(a), and 3(b)) and for case 3(c) where — — is actually consid- 

8o- k 

ered. Before this computation is made, however, the zeros in Jt 0 ,tfj of 
Pq[^(^,t),i//(5,t),a,t] (q = 1, 2, . . , l) need to be discussed. Given x(5,t) and ^p(a, t), 

let Pq(«,t) = p|$(«,t),^/(a,t),5,tj, and given a, let p q (t) = p q (a,t). As p q (t) passes to 
zero, sgn p q (t) may not necessarily change value. If r is a time point where 

p q (t) = 0, fictitious points occur when the first nonvanishing derivative of p q (t) from the 
left at r is of even order. These points occur where j5 q (t) has a local maximum or 

minimum (ref. 9). At such points, sgn p q (r~) = sgn p q (r + ). However, sgn p q (t) does 

change value for p q (r) = 0, and the first nonvanishing derivative from the left at t is 
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of odd order. These values of r, denoted by t*, are referred to as switching points of 
the switching function Pq(t) to distinguish them from zeros of Pq(t) where sgn Pq(t) 
does not change sign. Let the set of switching points of pq(t) be denoted by Sq(t*) and 
S(t*) denote the set of all switching points of all Pq(t) (q = 1, 2, . . . 1). Because the 
total number of switching points is assumed to be finite, the elements of S(t*) can be 
ordered such that 

S(t*) = (t* t* . . . t*) 


where t? + j > t? for i = 1, 2, . . . v. 

Cases 1, 2, 3(a), 3(b) .- In order to compute |=(a,tf) and when tj is 

fixed, p^a, tf) is derived by assuming that is known. The equation for 

then follows from observation of the equation for -p^a,tf). 


9a' 


When x(a,t), i^(a,t), and a are given, the function f is a piecewise continuous 
function of t with points of discontinuity occurring only on S(t*). Thus, x(a,t) is con- 
tinuous in t and satisfies the integral equation 


x(a,t) = x(a,t 0 ) + f f(x,u,a,s)ds 
x ' J t 0 


When x(a,t), \[/(a,t), and a are given, f (t) = f (x,u,a,t). Because f(t) is bounded on 
S(t*) 

x(a,tf) = x(a,t 0 ) + Jj( t ^f(x,u,a,s)ds 


= lim 
6-0 




f (s)ds 


where 

J(tf) = [t 0 ,tf] - S(t*) 

Q-P — — — — __ 

Let -g=(a,t) = 7 ~(x,u,a,t) for given x(a,t) and and a. Furthermore, because 

|^(a,t) is bounded and continuous on j(tf) 
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“ Sp.tf) + i 1 - e ) 

i=1 M*i*) 


~ dt? 

-f(t* + e)-i(5,t* + e ) 
8a 


+ C |L(a,s)ds 


l-l 


where the symbol means the sum over all Pq (q = 1, 2, . . . f) having tj” as 

«ltf) 


a common switching point. 


For arbitrary te£t 0 ,tf|, -g=(a,t) becomes the solution of the integral equation 

i\ i. fenSe-tf) 


i'l [_Pq(tj) 



where 


and 


fo 

(t < t*) 

H (‘ - ‘i) - \l 

\ J / 

(‘ > ‘f ) 


J(t) = [t 0> tj - S(t*) 


For a switching point t* of /5q(a,t) and a parameter vector a, an infinitesimal 
change da in a causes a change dpq(a,t*) away from zero. If pq(a + da,t* + dt*) 
is also to be zero, then 


dt 1 


= 1 


dp 

m ^(a,t*)da k 


k=l 


d Pq 

dt 


(a,t*) 
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where 

9 P, 


8 ak < “’ t * ) ' 


i=l 


+ 63^’“’^ 


and 


d 

dt 


-p q (a,t*) = ~ P q [^(a,t),j//(a, t ),a!,t| 


implies that 


The condition that the zeros of pq(a,t) in t change continuously with £% 

9 

k 




p q Ca,t) 


(q = 1, 2 , 1] k = 1, 2, . . . m) are continuous functions 

8 Po 


of t at t* for given a. Thus, in order to preserve this continuity, -^^(a,t*) must 


9 Pq 


vanish at any switching point where Pq(a,t*) = 0. In general, 


9a k 


<5,t*) 


replaced by 


d^ 1 


dtP- 1 



d P a 

dr< 5 -‘*> 


can be 


d p P r 


where p is the order of the first nonvanishing derivative 


dt? 




from the left of Pq(t) at t*. Because only switching points of Pq(t) are used, p can 
be considered odd. The change in a switching point t* with respect to a change in a 
parameter a k for a switching function Pq is then given by 


d 


P-1 


9t *-(a,t*) = - dtP 


9 Pr 


da k 




9 “k 


d p P, 


(14) 




dt* 




which, by assumption, is continuous at t*. By using equation (14) and by recognizing 


that sgn p n (t) is fixed over the intervals between the switching points, -E(a,t) can be 
11 da 

written as 
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(15) 


£1 

I [% + )-%r 

J'MnJtfl 




dt P-l 8a 


dP ^q/- 


>4 - *j*) 


+ C (?L + K | af 8u 8^ | 9f au | 9f 

Jj(t) \9x 8u dxjda 9u dip Sa 9u Sa Sa 


ds 


for given values of x(a,t), \f/( a,t ), -^(a,t), and a. 

da 

Computationally, this integral equation can be solved in the following manner. 
Given x(a,t), i//(a,t), ~(a, t), and a , integrate the differential equation 

' OQ! 

d/9x\_ / 9f | 9f 9u\9x | 8f 9u dip 9f_ 9u 9f 

dt \ da) \ 9x 9u dx) da 9u Sxp 9o! 9u da Sa 


(16) 


With l|( to ) = flM from to to tl- Cal1 thi s ||(«,ti ). Replace ||(«,ti") by 


#- i D + i |'( t r)- f (*n 

V5) 


dt p-l 9« { a ’h) 


d Pp 




dt 1 


and use this as the initial condition for the integration of equation (16) from tj + to tg - . 

Repeat the process until the desired t is obtained. The quantity f (t? + ^ - f |t? _ j is 
simply the "jump" that f (t) takes in going from tj*“ to t- c+ . 


Because g is subject to the same restrictions as f, 2^(a,t) satisfies the inte- 
gral equation 


n dt p_1 da ' 


= f§(a,to) + j [4n - s(*r| dPa Gt .i 


r(* • *1) 


f /M + -5s i^\9x + M M + 

Jj(t) \9x 9u dxjda \9x du dip) Sa 9u 9a da 


ds 


(17) 
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Computationally, equation (17) can be treated the same as equation (15). Given a 
and the solution of equation (2), equations (15) and (17) yield a simultaneous set of matrix 

integral equations for the determination of and -^(a,tf) for use in cases 1, 2, 

3(a), and 3(b). It can be noted that -^(a,t) and are piecewise continuous func- 

dot da 

tions over [t 0 ,tfj with discontinuities occurring on S(t*). 


9i h . 9Xf , x 

The initial conditions g^-(a,t 0 J and ■£^-(a',t 0 J (i = 1, 2, . . . n; 

k = 1, 2, . . . m) are to be determined from the nature of a in a particular problem; 
for example, if = !/'i(t 0 ) 


(t)- 1 (1=1) 

8 “l' 0 )'\o (1.2, ...n) 


dx 


and 837( l o) = ° 0-1.2. 


n; 


j = 1, 2, . . . n). 


Case 3(c) .- Let |g(a,t f ) and ||(a,t f ), given by equations (15) and (17), be denoted 

by X(5,tf) and ^(5,tf), respectively. In case 3(c), an extra term must be added to 
equations (15) and (17) so that they become 


and 




= *(“>•() + ^>*f) + 


respectively. The variables and J^a, tf^ are computed by treating tj as a 

parameter and by using equations (15) and (17). In order to solve for J=(a,tf) and 
the linear system 
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and 




must be solved. 

This system requires the inversion of the (2n x 2n) matrix 



Difficulty may arise when equation (18) is singular or if too many significant figures are 
lost during the inversion process. 


Critical Review 

The procedure presented provides several important advantages. Only a forward 
integration and a single matrix inversion must be performed to compute the correction 
vector given by equation (10). The matrix to be inverted is guaranteed to be nonsingular. 
The direction of the correction vector lies between the direction given by the gradient 
and the Newton- Raphson processes. An additional advantage is that a final-time correc- 
tion can be made an integral part of the process provided that the final time is an unknown 
parameter; that is, corrections in the final time can be computed at each iteration as a 
component of the parameter correction vector. Also, integral equations are available for 
influence matrices that describe the effect of a change in the parameters on the terminal 
conditions. 

The process also has some disadvantages. The fact that complete convergence can 
be obtained from an arbitrary choice of assumed parameters is not established. The 
technique proposed is basically a boundary -condition iteration scheme. Such schemes 
generally have sensitivity and convergence problems (ref. 10). The process does not 
generally eliminate such difficulties. Finally, the technique cannot be used in problems 
in which a singular control (ref. 11) might occur unless such an occurrence can be 
predicted. 
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EXAMPLE CALCULATION 


Problem Statement 

In order to exemplify the usefulness of the foregoing method, the two- point boundary - 
value problem arising from the application of the Pontryagin maximum principle (ref, 4) 
to determine fuel-optimal lunar- rendezvous trajectories is considered. The dynamic 
equations for rendezvous with a target in a circular orbit are developed in appendix B. 

The maximum principle is applied in appendix C. The result is a two-point boundary 
problem of the type considered. 


Application of Algorithm 

The foregoing policy of exhibiting all variables of a function is not continued 
because subsequent equations are quite lengthy and involved. For example, instead of 
writing p(x,if/,a, t), p or p(t) is written with p(x, i//,a,t) implied by the defining 
equation. 

From appendix C, the system of equations corresponding to equation (2) is 
f l = k l = x 2 


f 2 = x 2 


f 3 = k 3 

f 4 = x 4 
f 5 = x 5 



+ sgn p] t^ 2 ^ 2 R s 3 (x 1 + R Sx ) ^ 2r 


2iipxrj 


x 4 


y|l + sgn pj if/ 4 

J2 2 R s 3 (x 3 + R Sy ) 

2]jifx 

\£ 3 

X 6 


j/[l + sgn p] t// g 

^Rs ( x 5 + Rs z ) 

2\fpXrj 

^ JL \ 

\jx 3 

v\l + sgn p] 


2c 



S X 


+ fi 2 Rsv 


C0 2 xj + 


2ooxa 


- 2iox2 + w 
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where 


S2 2 R s V 2 3J2 2 R s 3 d 


\£ 3 

\£ 5 

+ 2 on // 4 


J2 2 R s 3 ^ 4 

3fi 2 R s 3 d 

j 

\S 3 

fr 5 1 

- 2 coi // 2 - 


« 2r s 3 * 6 

3fi 2 R s 3 d 


( X 1 + R s x ) - ^2 


>2,i 


0 0 V5E 

^6 = ^6 = '^5 

& 


\£ 5 


( x 5 + R s z ) 


g 7 = = | 


+ sgn p 




2 „ 2 
x 7 


d = ^2 ( X 1 + R s x ) + ^4 ( x 3 + R s y ) + ^ 6 ( x 5 + Rs z ) 


# = j|^ 2 2 + ^4 2 + 


i lx = Vfx, 


c 2 + R s x ) 2 + ( X 3 + R Sy) + ( x 5 + R S Z )' 


and 


R s (t) = yR Sx 2 (t) R Sy 2 (t) + Rs z 2 (t) (as given in appendix B) 


(*7 - *o) 


x 7 
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Because j/x is bounded away from zero, f = col(fj, ... £7) and 
g = col^gj, . . , grjj are continuous in x, xp, and sgn p and piecewise continuous 
in t with points of discontinuity occurring at the switching points of the switching func- 
tion p. 


In the notation of the algorithm, the boundary-value problem presented in appen- 
dix C becomes 


1 = ^i(to) 

e l = x i(tf) 

2 = ^2(M 

e 2 = *2 ft) 


"6 = M to ) e 6 = x e( tf ) 

with ^7(t 0 ) - i/'o normalized and t f e[t;p(t) = 0,p(t) S o]. Basically, a problem such as 
case 3 exists and is solved by the method of case 3(a). 

The equations corresponding to equations (15) and (17) for this problem are 


8x, 9xj pi 8xj. 1 9x4 . , 

8^ (t) = 55Jpo) + J to “taj" ds ’ ssyfo) = 0 a = 1, 2, . . . 6; i = 1, 3, 5). The jumps 

of f j (i = 2, 4, 6, 7) and grj as t passes through a switching point of p are 

y^t*) 


-sgn p(t*“) 


f— - (i = 2, 4, 6), - £ sgn p(t*-) (i = 7), and - - gn P ^* 1 


x 7 2(t*) 

for grj. Therefore, the remaining equations corresponding to equations (15) and (17) for 
j = 1, 2, ... 6 are 


9x 2 (t) 9x 2 (t Q ) 


90! i 


9d!i 


jp- 1 

sgnp(tig-)» 2 (tj;) d tP-l to j 


k=l 


3 


H 


dp^ 

dtP 


K) 




||(1 + sgn p) j 


to 


i/'s 


dip* 8^g 

dip 9 ^ 2^4 8^7 * 2^6 ^7 




x 7 # 3 / 


9(3! i 


1 


X 7 #' 




(Equation continued on next page) 
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3£2^R s ^xj + Rs x ] 

| 2 _ f2 2 R s 3 2 

^1 + 

— 

3J2 2 R s 3 1 

( X 1 + R s x ) 

( x 3 + R =y) 

\£ 5 

O + w 

vs 3 J 

8af j 

i x 5 J 


0X3 

doi\ 


8x 4 

+ 2u> - — + 
8orj 


3S! 2r 3 k * R %)( X 5 * R sz) 

S _ R 


tfx' 


8x 5 2 (1 + sgn P ^2 9x 7 

9of j x 7 2 ^// 8<a! i 


ids 


'8x, 


= 0 


9 x 4 (t) ^ 8x 4 (t 0 ) 


8a j 9a ] 


r 


d p-l 8p(t*) 

sgnp(^> 4 (tk)iFT-9^— 


H t 


( l - *J) 


pt r 

J|(l + sgn p)| 


2^4 8 ^2 t / 1 ^ 4 2 \ 8 ^ 4 9 ^6 

x ? \£ 3 8a j lx 7 \/i// x^sfip^J da j x 7 # 3 9a j 


2 3 ( x i + Rs x )(x 3 + H Sy ) 


3S2 R s 




8xi 3 x 9 

— i _ 2 cj — - + 

8(3! j 8 q! j 


3n 2 R S 3 (* 3 + _ R 5 s y) 




°V 2 

+ Ci>" 


8X3 

8q! i 


\Jx z 

J(1 + Sgn p)^ 4 8 x 7 
x 2 # ^ 


2 3 ( x 3 + R Sy)( x 5 + Rs z) 


3JTR S 


ids 


V* 1 


^5 

8Q!i 


(ta*(*°) - °) 


ap- 1 8 K*k) 

8x 6 (*) 8x 6 (t„) \ sgnp(t*-)^ 6 (t*) dt p-l 8 tt] 


r t 


0q? 3 9a j 


H(t 




k=l 




( t -‘iO 


s |(1 + Sgn p) 


(Equation continued on next page) 
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X 


'-* 2*6 ®*2 * 4*6 3 *4 I 1 _ *aiV* 6 


c 7 # 3 9o! i x 7 # 3 9 “3 \x 7 # x 7 # 3 y 9Q! 3j 


L 


2 n 3 ( X 1 + Rs x )( X 5 + R s Z ) 
\^ 5 


30 Z R S 


9Xj 

9of4 


3fi 2 Rs 3 ( x 3 + Rs y)( X 5 + R s Z ) 
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continuous at t,*. Time derivatives of are no t continuous at t,* because 

K 9Q-i k 

j J 

^ H(t - t*) is not bounded. In order to apply the algorithm, only the cases where the 

switching points of p(t) are simple zeros can be considered. For all such cases con- 
sidered, the number of zeros of p(t) were finite. 

Finally, because 'a - (a . . . . a^j and e(a,tf) = ^e^, . . . eg)', the matrix 
■^(3,tf) is the (6 x 6) array (dx^da^ (i = 1, 2, ... 6; j = 1, 2, ... 6). The mea- 
sure of terminal error is E(<5,tf) = 

i=l 



Results 

The algorithm and system equations were programed for the IBM 7094 electronic 
data processing system by using the Fortran IV language. Copies of the program are 
available on request from the Trajectory Applications Section, Langley Research Center, 
for the problem "Fuel Optimal Rendezvous" (program no. E1257). Integration was per- 
formed with fixed- step sizes of either 2 or 4 seconds by using a method with a fourth- 
order Adams-Bashforth predictor formula and a fourth- order Adams- Moulton corrector 
formula. 

The program was such that fixed-final-time and free-final-time solutions could be 
obtained. The program had the option of iteration at a fixed time or at a time when, after 
a specified number of coast periods have taken place, p = 0 and p g 0. The approach 
taken in constructing free-time solutions was to begin with a nominal and to compute suc- 
cessively fixed-time solutions for increasing values of the final time until, at such a time, 
a zero of p(t) was observed, which satisfied a prescribed number of coasts with pSO. 
This solution was then used as a nominal with p(t) = 0 and p ^ 0 as a stopping 
condition. 

The satellite orbital plane was placed in the xy-plane of the rotating system. (See 
fig. B-l of appendix B.) Examples were computed both with the vehicle launched from 
absolute rest from the surface of the moon in the satellite plane and from absolute rest 
from out of plane. Table I shows the values of the fixed parameters for classes of both 
examples. 
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TABLE I.- FIXED PARAMETERS 


Initial time, t 0 

Upper bound on thrust magnitude, y 

Initial mass, m Q 

Effective exhaust velocity, c 

Radius of moon, R v (to) 

Radius of satellite orbit, * R s 

Gravitational parameter of moon, p, 

Angular velocity of moon about axis of rotation, 


0 sec 

3504 lbm (1589.4 kg) 

285.5 slugs (4166.3 kg) 

9853.2 ft/sec (279.8 m/sec) 

5.707 X 106 ft (1739.4 km) 

6.1934 X 10 6 ft (1887.7 km) 

1.727 x 10 14 ft 3 /sec 2 

(48.9 x 10 11 m S/secty 
cd 2.66 x 10~6 rad/sec 


^Satellite is in an 80 n. mi. circular orbit. 


It was found that a workable set of b^ (i = 1, 2, . . . 6) and | X| for convergence 
was b^ = bg = b^ - 1, b2 = b^ = bg = 10, and |x| = 10. Except as noted, these values 
were used throughout. An increase in the value of |x| yielded slower convergence; 
whereas, a decrease was apt to produce divergence. Greater influence could be applied 
to the correction of the error e^ by increasing a particular bp that is, this error 
would be corrected more quickly than before, probably at the expense of the other errors. 
A similar statement could be made for the lack of influence on correcting e^ observed 
by decreasing bp 

In-plane results are presented in table II. When the vehicle was launched such that 
the satellite lead angle cp Q - cp v ° was 13.7° ( cp Q = 89°), the set of values 

^j(to) = 2-0 
^ 2 (t 0 ) = 3000.0 
<// 3 (t 0 ) = 10.0 
^ 4 (t 0 ) = 3000.0 
^tjO-o) = 0 
(to) = 0 

^ 7 ( t 0 ) = -3.3 x 10 6 

was found to yield a trajectory which continuously gained altitude with p(t) > 0 through- 
out; that is, no switching points were encountered. Then, \p 0 was reset such that p(t) 
went through a zero; that is, the vehicle began to coast at about 250 seconds. These 
values then produced a nominal set of values for - a g to which successive 
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TABLE II.- FUEL OPTIMAL IN- PLANE RESULTS 



Final time, 
tf, 
sec 

First coast 
time, 
sec 

Final burn 
time, 
sec 

Percent of 
initial mass 
at tf 


Unknown 

parameters 


10- 6 


^2 (W 

^3 (to) 

^4 (to) 

^ 7 ( t 0 ) x 10-5 

620.0 

313.5 

539.1 

51.07 

-0.19525 

2.43870 

3708.9 

10.6080 

3195.4 

1.12350 

720.0 

321.6 

657.9 

52.12 

-.18058 

1.54280 

3300.4 

7.5372 

2332.2 

-.97578 

850.0 

328.7 

800.0 

52.83 

-.17123 

.86229 

2991.7 

5.7089 

1801.4 

-.88233 

1000.0 

334.6 

959.2 

53.24 

-.16536 

.37087 

2766.5 

4.6379 

1477.4 

-.82356 

1150.0 

339.0 

1115.1 

53.43 

-.16189 

.05115 

2617.9 

4.0602 

1291.6 

-.78890 

1300.0 

342.5 

1269.1 

53.49 

-.15974 

-.16030 

2518.1 

3.7377 

1177.5 

-.76740 

1350.0 j 

343.5 

1320.8 

53.50 

-.15946 

-.18535 

2506.1 

3.6927 

1156.3 

-.76464 

a 1390.6 

344.3 

1361.6 

53.51 

-.15912 

-.21928 

2489.8 

3.6482 

1137.3 

1 -.76122 

l 


a Free-final-time solution, |p(1390.6)| = 0.432 X 10"®. 


application of the correction equation (10) at a fixed-final time of 620.0 seconds yielded 
the first entry in table II. Other solutions were computed by using the nearest obtained 
solution in final time as a nominal. At the 1350.0-second case, a zero of p(t) which had 
the property p(t) < 0 after one coasting period was predicted. With p(tf) = 0 and 
p(tf) < 0 as stopping conditions, the 1350.0-second case as a nominal, and |x| = 10 4 , the 
1390.6-second case resulted. 

For each of the solutions below the 1150.0-second entry, the trajectories remained 
below the altitude of the satellite orbit. Near the 1150.0-second final time, the tra- 
jectories, during the coast phase, began to reach an altitude which was higher than that of 
the target orbit. Such a property is referred to as "overshoot." Trajectories with and 
without overshoot are contrasted in figures 1 to 3. The trajectory without overshoot is 
the 850.0-second entry, and the one with overshoot is the free-final-time 1390.6-second 
solution. Arrows along the solid-line burning portions of the trajectories indicate the 
direction of the thrust vector. Dotted-line portions indicate coasting periods. The 
xy-axis system is represented as inertial because the total time of flight is such that the 
angular displacement of the moon about its own axis is less than 0.3°. 

Table in shows a sample iteration for the fixed-final-time 1300.0-second solution 
when the 1000.0-second solution is used as a nominal. Results for the vehicle launched 5° 
out of the target orbital plane are presented in table IV. Overshoot begins near the 
1200.0-second entry. In-plane results are used for beginning nominals. 
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Figure 1.- In-plane trajectories for <p 0 = 93.7° to final 
times of 850.0 and 1390.6 seconds. 


Figure 2.- In-plane trajectory for (p 0 = 93.7° to final time of 
850.0 seconds. 


1740 1790 1840 1890 km 

I F I I 

57 5.8 5.9 6.0 6.1 62x10° ft 

=□] — i — "r — i — fy 


* v °=807 


t7 = 344.3sec 


t#= 1390.6 sec ": 


Thrust portions of trajectory 
Coast portion of trajectory 

Direction of thrust vector 


t£= 1361.6 sec 


Figure 3.- In-plane trajectory for <p 0 = 93.7° to final time of 
1390.6 seconds. 


TABLE HI.- FUEL-OPTIMAL ITERATION SEQUENCE 
[<P 0 • <Pv 0 = 13 - 7 °; 0v° = °°; tf = 1300.0 sec] 


Iteration 

First coast 
time, 
sec 

Final burn 
time, 
sec 

Unknown parameters 

Error criterion 

to) 

t /' 2 ( t o) 

<P 3 (io) 

"o 

1 

E(or,t f ) 

*0 

334.6 

959.2 

0.370870 

2766.5 

4.6379 

1477.5 

0.101 X 10*3 

1 

360.8 

1296.1 

.401920 

2710.0 

3.9214 

1091.3 

.129 x 10 12 

2 

343.9 

1268.7 

.291960 

2734.2 

4.0687 

1256.0 

.776 X 10 9 

3 

342.3 

1259.2 

.145090 

2663.0 

3.9937 

1244.8 

.306 X 10 8 

4 

342.6 

1262.5 

.036295 

2611.7 

3.9098 

1221.7 

.228 X 10 8 

5 

342.5 

1264.7 

-.035214 

2577.6 

3.8512 

1206.4 

.726 x 10 7 

6 

342.5 

1266.5 

-.081362 

2555.7 

3.8110 

1196.1 

.234 X 10 7 

7 

342.5 

1267.7 

-.110730 

2541.7 

3.7844 

1189.3 

.596 x 10 6 

8 

342.5 

1268.4 

-.129310 

2532.8 

3.7672 

1184.9 

.137 X 10 6 

9 

342.5 

1268.8 

-.141050 

2527.2 

3.7561 

1182.1 

.290 X 10 5 

10 

342.5 

1269.0 

-.148460 

2523.7 

3.7491 

1180.4 

.563 X 10 4 

11 

342.5 

1269.1 

-.153140 

2321.5 

3.7446 

1179.2 

.105 x 10 4 

12 

342.5 

1269.1 

-.156100 

2520.0 

3.7418 

1178.5 

.204 X 10 3 

13 

342.5 

1269.1 

-.157940 

2519.2 

3.7400 

1178.0 

.338 X 10 2 

14 

342.5 

1269.1 

-.159100 

2518.6 

3.7389 

1177.8 

.120 X 10 2 

15 

342.5 

1269.1 

-.159840 

2518.3 

3.7382 

1177.6 

.167 X 10 

b 16 

342.5 

1269.1 

-.160300 

2518.1 

3.7377 

1177.5 

.600 


a Nomlnal tf = 1000.0 sec; x^t^) = 0.7616 X 10 6 , x 2 (t f ) = 0.52222 X 10 4 , 
x 3 (tf) = -0.12316 X 10 7 , and x 4 (tf) = -0.10558 X 10 5 . 

b Xi(tf) = -0.1824 !, x 2 (tf) = -0.23426, x 3 (t f ) = -0.36939, and x 4 (tf) = -0.21920. 


TABLE IV.- FUEL-OPTIMAL OUT-OF-PLANE RESULTS 
[to 0 = 93.7°; (p v ° = 80°; <V = 5°] 


Final time, 

First coast 

Final burn 

Percent of 




Unknown parameters 



sec 

time, 

sec 

time, 

sec 

initial mass 
at tf 

Xp 0 x 10-6 

^(to) 

*^2 0 "°) 

* 3 (to) 

^4 (to) 

>^ 5 (to) 

<^6 (to) 

^ 7 (t 0 )X 10 - 5 

620.0 

324.6 

528.4 

48.16 

-0.23492 

4.320300 

4718.2 

12.3090 

3937.9 

-6.6467 

-2179.00 

-1.51920 

720.0 

327.5 

645.1 

49.88 

-.20325 

2.805700 

3928.8 

8.1004 

2657.5 

-4.1913 

-1372.40 

-1.20250 

850.0 

331.7 

789.5 

51.13 

-.18462 

1.649800 

3368.9 

5.7369 

1921.5 

-2.8230 

-892.25 

-1.01620 

1000.0 

336.2 

945.0 

51.88 

-.17456 

.939600 

3032.4 

4.5158 

1525.5 

-2.0924 

-630.00 

-.91560 

1100.0 

338.8 

1054.8 

52.17 

-.17063 

.636920 

2889.6 

4.0662 

1371.4 

-1.8051 

-522.38 

-.87632 

1200.0 

341.0 

1158.6 

52.36 

-.16790 

.415080 

2784.7 

3.7718 

1246.0 

-1.6015 

-443.16 

-.83442 

1300.0 

343.1 

1261.6 

52'. 47 

-.16557 

.250510 

2706.5 

3.5776 

1186.9 

-1.3925 

-382.35 

-.82973 

1400.0 

344.9 

1363.9 

52.54 

-.16460 

.128460 

2648.0 

3.4511 

1130.3 

-1.3410 

-334.20 

-.81598 

1500.0 

346.5 

1465.8 

52.57 

-.16362 

.038713 

2604.4 

3.3715 

1088.1 

-1.2578 

-295.18 

-.80692 

1600.0 

348.0 

1567.4 

52.58 

-.16342 

.035952 

2601.4 

3.3406 

1062.6 

-1.2457 

-273.17 

-.80483 

a 1614.3 

348.3 

1582.0 

52.58 

-.16336 

.032755 

2599.6 

3.3362 

1058.9 

-1.2425 

-269.84 

-.80429 


a Free-final-time solution, |p(1614.3)| = 0.461 X 10"6. 
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Running time for all cases on the IBM 7094 electronic data processing system was 
approximately 7 minutes. In programing this example, the primary objective was to 
decide whether the method could be applied to such a nonlinear two-point boundary-value 
problem and not necessarily to write a program giving solutions in a minimum of com- 
puter time. Time-consuming subroutines were included to test for conditions leading to 
numerical instability (overflow, underflow, and so forth) and to determine the nature of 
the zeros of the switching function. These subroutines and the large number of equa- 
tions involved account for the rather long computer time. No sensitivity or convergence 
problems were found in any of the cases considered. 

CONCLUDING REMARKS 

A successive approximation procedure for attacking a class of two-point boundary- 
value problems which frequently occurs in indirect optimization theory has been pre- 
sented. Basically, the boundary- value problem was one in which the optimal- control law 
was piecewise continuous and in which there were a number of system parameters to be 
determined to meet an equal number of terminal conditions. An iterative logic was 
developed in which an assumed set of parameters would be improved upon so that, by 
repetitive use of a correction formula, a monotonic decreasing sequence of values of a 
positive definite function that measures the terminal errors was produced. 

The procedure provided several important advantages. A forward integration and a 
single matrix inversion must be performed to compute the correction vector. The matrix 
to be inverted was guaranteed to be nonsingular. The direction of the correction vector 
was found to lie between the direction given by the gradient and the Newton- Raphson pro- 
cedures. An additional advantage was that a final-time correction could be made an 
integral part of the process provided that the final time was an unknown parameter; that 
is, corrections in the final time would be computed at each iteration as a component of 
the parameter correction vector. Integral equations were derived for influence matrices 
that describe the effect of a change in the parameters on the terminal conditions. 

The process also had some disadvantages. The fact that complete convergence 
could be obtained from an arbitrary choice of assumed parameters was not established. 
The technique proposed is basically a boundary- condition iteration scheme. Such schemes 
generally have sensitivity and convergence problems. The process does not generally 
eliminate such difficulties, but none were found in the example considered. Finally, the 
technique cannot be used in problems in which a singular control might occur unless such 
an occurrence can be predicted. 

In order to demonstrate the usefulness of the procedure, solutions were obtained to 
the two-point boundary -value problem resulting from an application of the Pontryagin 
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maximum principle to obtain fuel-optimal lunar- rendezvous trajectories for a target in 
a circular orbit. Fixed- and free-final -time solutions were computed for planar and 
nonplanar situations. Running times on the IBM 7094 electronic data processing system 
were on the order of 7 minutes. In programing this example, the primary objective was 
to decide whether the method could be applied to such a nonlinear two-point boundary- 
value problem and not necessarily to write a program giving solutions in a minimum of 
computer time. Time-consuming subroutines were included to test for conditions leading 
to numerical instability (overflow, underflow, and so forth) and to determine the nature of 
the zeros of the switching function. These subroutines and the large number of equations 
involved accounted for the rather long computer time. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., May 8, 1968, 

125-19-04-01-23. 
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APPENDIX A 

PROOF OF LEMMAS USED TO ESTABLISH THEOREM 1 
The lemmas used to establish theorem 1 are now proved. 

Lemma 1 

Unique solutions 6a°(v) and X(v) <0 of the system 

— nl|2 ~ ^ 


||6q'°|| = 

6 "° = "[^(" 0>tf ) B H^° >tf ) " XI 


^(a°tf)Bi(a°,tf) 


(Al) 


exist if v is sufficiently small. 


Proof: Note that tf^B is a real symmetric matrix and can there- 


da' 


da 


fore be diagonalized (ref. 12). There exists an orthogonal matrix A,A T = A \ which 


operates on -^(a°,tf)]3 -^=-(a°,tf) to yield 


A ~|(a°,tf)A' = diag(Xi) 


(i = 1, 2, . . . m) 


where the are the eigenvalues of -^(o! 0 ,tf)B |=(a! 0 ,tf). If rj is an arbitrary 


m-dimensional column vector, then 


9cr 


da'' 


where \/b = diag^bj because B = diag(bj) (i = 1, 2, . . . m). Thus, 

■^|r(a' 0 ,tf)B |=(a°,tf) is nonnegative definite and, therefore, has nonnegative eigenvalues 
(ref. 13); that is, Xi = 0 for all i = 1, 2, . . . m. Therefore, the inverse of 


^[a°, tf)B -^1(5°, tf) - XI exists because X < 0. When 


95' 


and 


90 ! 


c = -A^(a! 0 ,tf)Be(a 0 ,tf) = col(ci) 


(i = 1, 2, . . . m) 


6v = A6a° 
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this transformation reduces 

■^lr(a 0 ,tf)B -p(a°,t f ) - XI 
9a' ’ v 9a v v 

through 


A'6v = - A ~(a°,tf)Be(a°,t f ) 
to 



6 a° = - — (a °, tf) B e (a °, tf) 
9a ' 


diag(Xi - X)6v = c 


Because X < 0 


and 


6v = diag c 

Xi ~ X 


m 2 

||6a°||^ = 6a° • 8a° = 6v • 6v = V — — — — = 

i=l ( X i " X ) 


or 


m 9 

y c i 
i=l ( x i + l X D 2 


Assume that Xj # 0 for all i = 1, 2, . . . m. If 

m • >2 


^2 > V / c iV 

" = i 4 ^ 


(A2) 


no real negative X exists because the expression 


is strictly decreasing for increasing 



i + I x l) 

For 


2 


v 


2 



a unique X which satisfies equation (A2) exists. If some Cj vanish, these terms in 
equation (A2) vanish independently of Xi and the same arguments hold for the reduced 
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equation. If, for i » j (j = 1, 2, . . . m), Xj = 0 but Cj 4 0, then equation (A2) 
becomes 

c.2 _ m - c .2 

and solutions in X exist for all v. Finally, if all Cj vanish, a solution exists only for 
v = 0. 


Lemma 2 

The solutions 6 a°(v) and \(v) <0 of the system (eq. (Al)) maximize the absolute 
value of 

SE(a°,t f ) = 6 Be(3°,tf) • 63° + | 63° • -^lr(a 0 ,tf)B ||(a°,tf)6a° (A3) 

subject to the conditions || 6«°|| = v 2 and AE(3°,tf) < 0. 

Proof: Note that the inequality condition on ||63°|| can be replaced by an equality 

2 

condition through the introduction of a real variable /3 because ||63°|] = v 2 implies 

and is implied by the existence of a jS such that 63° • 5 3° - v 2 + / 3 2 = 0. Then, 
|AE(3°,tf)| must be maximized with respect to the choice of 63° and /3 subject to 

Condition (a): 

||63°|| 2 - v 2 + /3 2 = 0 

and 

Condition (b): 

SE(a°,tf) < 0 

A Lagrange multiplier X (ref. 9) is introduced, and 6a° and /3 are chosen such that 
the augmented relation 

|AE(a°,tf)| = |AE(3°,tf)| + |(||63°|| 2 - v 2 + / 3 2 ) 
is maximized. If condition (b) and equation (A3) are taken into account 

| AE(3°,t f )|* = - ^(3°,t f )Bi(3°,t f ) -630-1 5 -o . -^(3°,t f )B ||^o,t f )63 0 
+ ^(s3 0 * 63° - v 2 + /3 2 ) 


I 
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Necessary and sufficient conditions that 
to 6 a° and /3 are 

Condition (c): 


2E(a! 0 ,tf)| be maximized with respect 


8 SEfe°, t f ) 

J — = 0 


8 k 


and 


Condition (d): 


Sa ( 


AE(a°,t f ) 


8k 2 


is negative definite 


where k = . Condition (c) yields the vector equation 

\ P I 


{- ~(5°,tf)Be{5°,t f ) - ||(5 0 ,tf)B ||(5° t,)55° + X65°\ / o\ 

\ Xf ) \o) 

and condition (d) yields the matrix condition 

(- If 5 ”'!) + XI > 5 


8a v 




< 0 


o' 


, V 

A 


If A diagonalizes •^| r (a 0 ,tf')B |£(a°, tf), then G = diagonalizes 

da v J y \n 1/ 


a 2 SE(a°,tf) 


or 


AE(a°,t f ) 

r / 

_ G f = 

ail 2 

1 


/diag(x - \ij, O' 


Because G is nonsingular, examination of 


AE(a°,t f ) 


8k 2 


G’ 


for negative definiteness is equivalent to the examination of 


I 
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9 2 |AE(a°,t f )|* 

Because Xj. s 0 (i = 1, 2, . . . m), — 1 ^ is negative definite for arbitrary 

ak 2 

Xj if and only if X < 0. From X/3 = 0, /3 = 0, whereby 

5a° • 5a° - + /3^ = 0 


and 




5o»° : 0 


yield equation (Al). 


Lemma 3 

The quantity AE^a°,tf), given by equation (A3), is negative definite if 5a 0 for 
6a° satisfies equation (Al). 

Proof: From equation (Al) 

65° ■ |L|S°tf)B5(5 0 ; tf) = X(65° ■ 65°) - 65° • ^(5° t f )B ||(5°,t t )65° 
which upon substitution into equation (A3) yields 

AE(5°, t t ) = -|*|(6S° • 65°) - I 65° . |t(5°,t()B §(5°, t t )65° 

- f - 

Thus, 2E(a! 0 ,tf) is negative definite in 65° because |>|* 0 and -^•^°,tf)B |=(o! 0 ,tf) 
is nonnegative definite. Therefore, AE(a°,tf) = 0 if and only if Sa° = 0. 
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DYNAMIC EQUATIONS FOR LUNAR-RENDEZVOUS PROBLEM 

Dynamic equations are developed for a space vehicle which seeks to rendezvous 
with a station in a circular orbit in the vicinity of the moon. 

The vehicle is a one- stage rocket, treated as a point mass, with bounded thrust 
magnitude, and the controls are the magnitude and direction of the thrust vector. 

Let x, y, and z be Cartesian coor- 
z dinates of a rotating axis system located in 

the center of the moon with the z-axis through 
the axis of rotation of the moon. The geom- 
etry is represented in figure B-l. 

The vector R s (t) is from the origin 
to the space station. Let R v (t) be the 
y instantaneous vector from the origin to the 
vehicle and <x> be the angular velocity of 
the moon about its axis of rotation. By 
assuming that 

A m(t) = - £ (m(t 0 ) = m 0 ) 

Figure b-l- Rotating axis system. the dynamic equations of motion for the sta- 

tion and vehicle are 


(Bl) 


where 

CO = kco 

R s = i R s x + i R Sy + kR Sz 
Ry = i Ry x + j Rvy + k R y z 

i,j ,k unit vectors along x-, y-, and z-axes, respectively 
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m(t) total vehicle mass 

\i universal gravitational constant multiplied by mass of moon 

R v = (R v • Rv) 1 ^ 2 
R s = (R s • R s ) l/2 

T thrust control vector of vehicle 

T magnitude of T 

c effective exhaust velocity of vehicle rockets 


The thrust vector is related to the xyz-axis system by 


T = i^T cos e c cos cp + f(T cos 0 C sin <p c ) + k^T sin 0 C ) 


as shown in figure B-2. 

The vector R s (t) can be found at any instant in the rotating coordinate system by 
integrating its differential equation with the appropriate initial conditions or by a mapping 
process. Because the satellite is assumed to be in a circular orbit, it moves in its 
orbital plane at a constant distance R s from the center of the moon with a constant 
angular velocity £2 = . Consider an inertial XYZ-axis system fixed in the 

center of the moon such that, at the initial time t 0 , it is alined with the rotating 
xyz-axis system. In this framework, the station can be pictured as in figure B-3. 


The angles i Q and 0 o define the normal and line of nodes, respectively, of the 
target orbital plane relative to the inertial system. The x' and y' axes define the 
orbital plane of the target. If, at t Q , the target is in the position 
(x',y') *(k s cos <p 0 ,Rs sin and moves clockwise from R s , then 


RspWjy'^z’ft)] 



(B2) 


Therefore 

R S (X,Y,Z) = Tiiis(x’,y’,x’) 


(B3) 
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2 



Figure B-2.- Reference axis system for control 

vec tor. Figure B-3.- Station viewed in inertial axis system. 


where 



-cos i Q sin 0 O 
cos 6 Q cos 0 O 
sin i Q 


sin o 0 sin 9 0 
-sin i Q cos e Q 


cos L 


(B4) 
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or 


Also 


with 


and 


or 


R s (t) = R s < 


COs|w(t - t 0 ) - 0 qJ 
sin|o) (t - t 0 ) - 6>o] />cos|^ 0 - S2(t - t Q )j 


cos u Q sinjw(t - t Q ) - flo] 

+ /cos c 0 cos|w(t - t Q ) - tfjjl sinjV Q - fl(t - t 0 )jj 


sin i r 


= T 2 (t)R s (X,Y,Z) + T 2 (t)R s (X,Y, Z) 


T 2 (t) = co 


-sin co(t - t Q ) cos co(t - t 0 ) 0 

-cos co(t - t D ) -sin w(t - t Q ) 0 

0 0 0 


r sinj> 0 - 0(t - t 0 )j 
R S (X,Y,Z) = R s OT^-cos^ 0 - 0( t - t Q ^ 


dR s (t) 



-sin£co(t - t Q ) - ©oj^co + S2 cos (, Q ) 
cosjco(t - t 0 ) - 0<^o> + fi cos t 0 j\cosj^c> 0 - J2(t - t 0 Jj 


-J2 sin c r 


J 


(B6) 


(Equation continued on next page) 
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r x \ 

COS^w(t - t 0 ) - ejj (w COS 6 0 + O) \ 

+ < -sinjw(t - t 0 ) - 0 o ](o> cos l 0 + J2)\sinj^ 0 - S2(t - t 0 )j (B7) 


Thus, the position and rate of Rs(t) can be obtained by specifying R s , l 0 , 6 0 , 

and cp 0 at t 0 and by using equations (B6) and (B7). The initial value R v (to) can be 
specified by 

R v (t 0 ) = R v (t 0 )(i cos 9 V ° cos (p y ° + j cos 6 y ° sin cp y ° + k sin 0 y °j 

where 9 y ° and cp y ° are shown in figure B-4. Also 

2 Rv(to) = i Rv x ( t o) + J Rvy(^-o) + kRy z (to) (B8) 

v / Rewriting equation (Bl) in vector- matrix 

notation yields 



Rv(t 0 ) U 4 U R ^ 9 _ 

! Rv = M ' 2SRv " s ^ 



= R v °; R„(t 0 ) = Rv“) 


Rs “ "M — t “ 2 SR s - S Rg 
R s ^ 


Figure B-4.- Initial orientation of vehicle with 
respect to rotating axis system. 


(Rs(*o) = R s°; R s(to) = R s°) j 


where 


Ro = R. 
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" u l~ 


cos Be cos cp c 

u = 

u 2 

= 

cos Be sin cp c 


_ u s_ 


sin Be 


and 


S = 


u 4 = T 


0 -co 0 


co 0 

0 0 


In terms of the relative distance ? = R v - R s 


' - silij - R2f - n2 ( f + R s)( 


R, 


s -l)-2Sr-S 2 r 


r + 


R sf 


where 


A 

r + R s = (r + R s ) • (r + R s ) 


(BIO) 


Because the maximum principle is used for optimization purposes in appendix C, the 
state-vector notation is now employed. Let 


and 



x 6 = r z = x 5 
x 7 = m(t) 

x 8 = t 


With R s and Rg regarded as explicit functions of time by equations (B6) and (B7), 


write 
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. u^Mu 
v = — — 


u 4 

x 7 = -“ 


x 8 = 1 




+ Y (v,x 8 ) (v(t 0 ) = v Q ; v(tf) = 0) 

(x 7 (t Q ) = m(t 0 )) 'f 
( x 8 (to) = t Q ) 


(Bll) 


where 


f2^R 2 » 

y(v,x 8 ) = — -|Nv + MR s (Xg) + f2 2 MR s (Xg) + (n’ + 2coK + w 2 l)v 

||Av + R s || L J 

v = COl^Xj, . . . XgJ 


launch time 


tf final rendezvous time 



0 

0 

0 

1 

0 

0 


0 

0 

0 

0 

0 

1 


0 0 0 0 
10 0 0 
0 0 0 0 
0 0 10 
0 0 0 0 
0 0 0 0 


0 0 
0 0 
0 0 
0 0 
0 0 
1 0 



0 0 0 0 0 

0 0 10 0 
0 0 0 0 0 

-1 0 0 0 0 

00000 
0 0 0 0 0 
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0 

1 



0 

0 


0 0 
0 0 
0 0 
0 1 
0 0 
0 0 


0 0 
0 0 
0 0 
0 0 
0 0 
0 0 


0 

0 

0 

0 

0 

0 



0 0 0 0 0 
0 10 0 0 
0 0 0 1 0 


The act of rendezvous requires that the vehicle and station have the same position 
and velocity at tf; hence, the condition v(tj) = 0. In addition, u^ s y, where y is the 
largest value obtainable for the thrust magnitude. 
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NECESSARY CONDITIONS FOR FUEL-OPTIMAL RENDEZVOUS 


Given the system of equations (Bll), establish necessary conditions that the con- 
trol functions u and u 4 drive v(t) from v(t 0 ) to v(tf) = 0 in such a way as to 
A ftf u 4 

minimize x Q ^tf) = j — dt = m(t G ) - m^tf). These conditions readily follow from the 
*o 

Pontryagin maximum principle (ref. 4), 

Before the maximum principle can be stated for this particular problem, the vari- 
ables x Q , H, and i^ k (k = 0, 1, . . . 8) must be defined as follows: 


-'n 


X n = 


u 4 


(x 0 (to) = o) 


H 


O 

■= 2 V: 


} (Cl) 


k=0 


9H 

i!/, = - — 

k 9x v 




Through equations (Bll), the equations for (k = 0, 1, . . . 8) are 


where 


^o =° 


fi = - 

8Y \ ,x 8 ) , 

9v 

^ 7 = 

u 4 h • Mu 

x 2 
x 7 


8Y(v,x 8 ) • 

00 

II 

9xg 


h = col^p rp 2 , . . . 


The Pontryagin maximum principle can then be stated as follows: Let Uj 

3 

(i = 1, 2, . . . 4), where ^ u^b 1 and 0 % u 4 % y, be controls which transfer xj(t Q ) 

i=l 
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to (j = 0, 1, ... 8). In order that Uj minimize x 0 (tf), it is necessary that 

there exist a nonzero continuous vector with elements (j = 0, 1, . . . 8) as deter- 
mined by equation (Cl) such that: 

(1) For every t (t Q £ t g tf), the function g^Xj,^j,u^, for fixed Xj and i//j, 
attains its maximum at the point uj = uj(t) ; that is, Hp//j,Xj,u^(t)j = M^j,Xj 


If equations (Bll), (Cl), and H|^j,Xj,uj(t)j = M(^j,Xj) are satisfied, then \p Q £ 0 and 
M|xj(t),i//j(t)J = 0. 

(2) Because x^t^ and Xg(tf^ are unconstrained, the part of the maximum princi- 
ple known as the transversality condition requires that ^(tf) = ^g(tf) = 0. Substitution 
of xj (j = 0, 1, . . . 8) into H gives the equation 


H = u. 




u = 


If M'h * 0, then the u which maximizes g and satisfies u • u = 1 is 
^ ^ The assumption M'h = 0 over a finite interval in [t 0 ,tfj leads to the con- 


tradiction i//j^tf) = 0 (j = 0, 1, ... 8) and, therefore, cannot occur on an optimal tra- 
jectory. If M'h = 0 at isolated points of [t 0 ,tf|, then the continuity of i^j(t) implies 

that, if t' is such a point, u(t') = ,. M — Xr. Then, by using the optimal direction 

||M'h(t’-)|| 

M ^ g becomes u ^ - —^ ~~^~ + + ^ * Y (^’ x s) + ^8‘ The u 4 which maxi_ 


u = 


|M’h|’ ~ 


mizes g, if p(t) = 1 


^7 ^o r i 

— — + — vanishes only at isolated points within |t 0 ,tfl, is 


u 4 = 


Y 

0 


(p> 0) 
(P< 0) 


or 

u 4 = ||l + sgnp] (C2) 

Situations in which p = 0 over a finite interval in jjto>*f] are referred to as being 
singular. Such cases may always occur in a general problem when the control enters 
linearly in H. In this case, the coefficient of the control must be examined to determine 
whether there is an admissible control rendering it identically zero. Such a control is 
termed a singular control. The existence of a singular control brings into account the 
difficult question of uniqueness of optimal controls. Necessary conditions for singular 
controls to be optimal are given in reference 11. 
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In general, the existence and optimality of a singular control rendering p = 0 
within [t 0 ,tfj should be examined. However, because the purpose of the fuel-optimal- 
rendezvous problem is to exemplify the use of the algorithm in the solution of a class of 
two-point boundary- value problems, singular solutions are not considered. All solutions 
obtained are for the nonsingular control law 

u = u 4 u = £[l + sgn pi M jj.j (C3) 

4 2L J ||Mh|| 

Substitution of u into H yields 

= 2 C 1 + Sgn P ] P + ^ ' ^’ X 8) + ^8 

At tf, where v(tf) = 0 and i//g(tf j = 0, M(tf) = 0 yields 

|[l + sgn p]p(t f ) = 0 


Thus, at tf, p(tf) S 0. Because, theoretically, a coast into a rendezvous cannot be per- 
formed, p(tf) = 0 (p(tf) S 0) . This property classifies tf ; namely 

tfe[t;p(t) = 0,p(t) £ 0] 

The equations to be satisfied along an optimal trajectory take the form 


X ° = 2c\ } + Sgn P ] 

* = ID + sgn p ]iiS 


“\ 

(x 0 (to) = 0) 


MM'h 
M'h x- 


fi 2 R 3 


o[Nv + MRgfXg)! + S2 MR s (x 8 ) 

.. .. 7 j|Av + Rs|j L /J 

+ (n' + 2coK + cu^Lijv ^v(t 0 ) = v Q ; v(tf) = 0) 


= _ J_[l + Sgnp ] 


k Q = l 




( x 7 (to) = m 
(x 8 (t 0 ) = t Q j 

(^0^0) " °) 


f2 2 R s 3 N'h 30 R s { 

Nv + MRg^Xg^ 

• S J 

[A ', 

Av + R s 

3 

Av + R s 

5 1 


Av + R s 


N] 


(n + 2coK' + co 2 L’)h 


(h(t Q ) undetermined) J 
(Equations continued on next page) 


> (C4) 
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Because t Q is interpreted to be the launch time 

(to) = -jyp(to) + h(t 0 ) • y(v 0 ,t 0 ) 

implies that M(t 0 ) = 0 whereby M(t) s 0 on [to,tfj which guarantees that i//g(tf) = 0. 
Thus, i// 8 (t 0 ) can be determined to eliminate the condition i//g(tf) = 0. If Xg is 
replaced by t, then the variables x Q , Xg, and can be eliminated and computed, if 

desired, after v(t), x^(t), \p Q , ^(t), and h(t) have been found. 

Thus, a two-point boundary-value problem occurs which requires the determination 
of i/^j(t 0 ) (j = 0, 1, ... 7) such that at tf, Xj(tf) = 0 (j = 1, 2, ... 6) and ^(itf) = 0. 
A reduction in the number of parameters i//j(t 0 ) and terminal conditions can be made by 
observing that the combination i//^(t) - \p Q satisfies the same differential equation as 
t) and, therefore, can be computed collectively. If a value of \prj(t 0 ) - can be 
set for which i//j(t 0 ) (j = 1, 2, . . .6) can be found such that v(tf) = 0, but t^(tf) * 0 
necessarily, then new values of 1^/7 (t 0 ) and \J / 0 can be computed for which 1^7 (tf) = 0 
and the combination remains the same. These new values are found by replacing 1^7 (to) 

cl M'h(tf)|| 

by ^ 7 (t 0 ) - ^ 7 (tf) and xp Q by xp Q - 1/^7 (tf). From p(t f ) = 0, ^ Q = ~ x 7 (t f y ~ and 
verifies that £0. 

Finally, the problem is to find ii(to) such that v(tf) = 0 with 1^7 (t Q ) - normal- 
ized and tfe|t;p(t) = 0,p(t) £ oj. Such a solution yields a trajectory which satisfies the 
necessary conditions for free-final-time fuel optimality. 
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In conclusion, if tf is fixed, does not necessarily equal zero (ref. 4) in 

which case ^Pq(H) can be adjusted to satisfy M(tf) = 0 and, thus, eliminate the necessity 
of p^tf) = 0. If a solution can be obtained with p(tf) arbitrary but ip Q 5 0, then the tra- 
jectory satisfies the necessary conditions of the Pontryagin maximum principle for 
fixed-final-time fuel optimality. 
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